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SUMMARY 

This report extends the results of Vinti and Izsak and pre- 
sents a'computational procedure designed specifically for Izsak’ s 
second-order solution of Vinti’ s dynamical problem. With this 
procedure, the coordinates and velocity of an unretarded satellite 
can be obtained from a knowledge of its initial conditions r and v. 

In this procedure, the derivation is given for the complete 
set of six canonical constants from initial conditions. Three of 
these have been determined by Vinti and the remaining three by 
the author. All six of them are assumed known in Izsak’ s solution. 

This report includes an adaption of a Newton- Raphson itera- 
tion scheme specifically designed to solve a certain system of 
nonlinear equations introduced by Vinti for the purpose of numer- 
ically factoring a certain quartic equation. The solution by this 
method can be used instead of certain infinite series to obtain 
Izsak’s elements a and e. An example is included to illustrate 
how these elements may be obtained by the Newton- Raphson 
method. 

Appendix B gives the derivation of exact expressions for the 
components of velocity in Vinti’ s accurate intermediary satellite 
orbit using Izsak’s orbital elements. The derivation is one of the 
necessary steps in comparing such a method with others. 
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A SATELLITE ORBIT COMPUTATION PROGRAM 
FOR IZSAK’S SECOND-ORDER SOLUTION OF 
VINTI’S DYNAMICAL PROBLEM 

by 

Raymond V. Borchers 
Goddard Space Flight Center 


INTRODUCTION 

This report provides a computational procedure for determining the orbit of an arti- 
ficial satellite in the earth’s gravitational field. The procedure is based on Izsak’ s second- 
order solution of Vinti’s dynamical problem (Reference 1). This computing procedure 
differs from many other methods in that the potential function is included in an analytic 
solution of the equations of motion. This is advantageous because the difficulties associ- 
ated with the slow convergence or divergence of some series expansions used in orbit cal- 
culations are avoided; also the problem of small divisors is avoided. Another advantage 
of this procedure is that it does not involve several multiplications of Fourier series, a 
task common to certain satellite programs. Although Fourier series are well adapted to 
numerical computation, it is certainly desirable from the standpoint of machine storage 
and computing time to minimize the total number of such series. In many satellite theories, 
Fourier series are used from the very beginning to obtain successive approximations of 
different orders to the solution. The use of Vinti’s potential minimizes the use of pertur- 
bation theory; Izsak (Reference 1) states that the oblateness perturbations which are not 
accounted for by Vinti’s potential can be treated by a first-order method, that is, without 
multiplications of Fourier series. 

As was pointed out by Izsak (Reference 2) it is advantageous for several practical 
purposes to have satellite orbits with very small eccentricities. Since the eccentricity 
never appears as a divisor, this procedure is valid for arbitrarily small values of c or 
e = o. However, we must avoid polar orbits and orbits which have inclinations of less than 
2 degrees. 

Vinti (Reference 3) found an axially symmetric solution of Laplace’s equation in oblate 
spheroidal coordinates which may be used as the gravitational potential about an oblate 
planet. This potential, which leads to separability of the Hamilton- Jacobi equation, is a 
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remarkable approximation to the actual gravitational field of the earth in that it fits the 
zeroth and second zonal harmonics exactly and accounts for over half of the fourth zonal 
harmonic. Naturally, the oblateness perturbations are only a part of the factors which 
affect the satellite motion. Other perturbations not accounted for in this procedure are 
the effects of the odd harmonics, the residual fourth harmonic, the lunar-solar forces, 
and aerodynamic and electromagnetic drag. 


MATHEMATICAL PROBLEM 


In Hamiltonian form, the equations of motion of a dynamical system of n degrees of 
freedom assume the forms 


d Pj _ _ m 

dt " 3q ; 


dq i _ 

dt dP; 


1 


1,2, " ' , n 


> 


(i) 


where H(q,, q 2 , •••, q n ; p,, p 2 , •••, p n ; t) is the Hamiltonian function (in which time ap- 
pears explicitly) of the system with n generalized coordinates q,, q 2 , • • • , q n and the con- 
jugate momenta p,, p 2 , •••, p„. 


Solving the Hamilton-Jacobi equation 


<jw l aw aw aw \ _ 

at + H r 1 i’ q 2’ q n ; aq , • aq 2 ’ Jq n y 0 


where W is Hamilton's characteristic function, is equivalent to solving the Hamiltonian 
equations of motion (Equation 1). If it is possible to separate the variables in the Hamilton- 
Jacobi equation, then the solution can always be reduced to quadratures. 

Vinti's dynamical system belongs to a class of systems which are scleronomic, con- 
servative, and holonomic. Furthermore, it belongs to a class of dynamical systems which 
are said to be of Stackers type. The separability properties of the Hamilton-Jacobi equa- 
tion of the form solved by Vinti follow from certain conditions determined by Stackel. The 
separability of the variables occurs only in certain coordinate systems. 
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The oblate spheroidal coordinate system is related to the geocentric rectangular co- 
ordinate system by 


x - l/{p 2 + c 2 ) Vl w CT 2 COS a , 

y = ^ p 2 + c 2 fl - a 2 sin a, 

z - pa t 

r - }^p 2 + c 2 (l “ cr 2 j , 

where x, y, and z are rectangular coordinates; r is the geocentric distance of the satellite; 
p , cr, and a are the coordinates in the oblate spheroidal system; and c is a constant defined 
by Vinti T s expression 


c2 - J 2 a E 2 ' (2) 

where J 2 is the coefficient of the second-degree Legendre polynomial in the earth’s force 
function F. The quantity F is expressed as 


F = 



r 



n= 1 


where S is declination of the satellite, a E is the earth’s equatorial radius, and p is the 
product GM where G is the gravitational constant and m is the earth’s mass. 

The potential which Vinti obtained in oblate spheroidal coordinates is 

y = 

p 2 + C^cr 2 


Similarly, the Hamiltonian and Lagrangian are 

H = 4u 2 

1 P 

L = \ U 2 + - 
* P 


2 


PP ( 
+ c 2 cr 2 


±P 

+ C 2 <7 2 
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where the speed of the satellite u is found from 

p 2 + C 2 Cr 2 . p 2 + C 2 Cr 2 , 


u 2 = - 2 

p 2 + C 2 


+■ 

9 1 - ^ 


a 2 + (p 2 +c 2 ) (l “ a 2 j a 2 


The generalized momenta are defined by 


db 

‘ P dp 


db 

do 


db 

do 


The Hamiltonian does not contain the time explicitly, so the Hamilton- Jacobi equation 


is 


2 (p 2 + c 2 a 2 ) 


/ \ /aw\ 2 

('° 2 + c2 )(^) + l 1 


d W 






p 2 + c 2 a 2 


= fi 


where, in the limit as c 2 -Oof Keplerian motion, fi is the total energy in the orbit; fi is 
always negative. 

The implicit equations of motion for the determination of p, a, and a are (Reference 1) 


d W 


c9W 

dc 


CP p 2 dp „ fi 47 cr 2 da 

1 ^ + C 2 I -7== ‘ t - t , 


L = f p p 2 < ] p , 2 f 

r " J VPGT + c J, Yom 


**r 


da 


d P 

/Wp) T do Vot 


gw 2 ~ a dp 

<?G C G J^ (p2 + c 2 ) /Pip! 


5 i: 


da 


(l “ a 2 ) VQ(a) 


+ a 


fi , 



where 


5 


P{p) - 2h p 4 + 2 pp 3 ~~ (c 2 “ 2c 2 hj p 2 + 2c 2 pp ~ c 2 {c 2 ~ G 2 j , 

Q(a) = - 2c 2 ha 4 - (c 2 - 2c 2 \\j cr 2 + ^c 2 - G 2 ^ . 

The symbols fi, c, G, - t, and fi are a canonical set of constants of integration. In the 
limit as c 2 -0 of Keplerian motion, the canonical constants have the following meanings: 

fi total energy in the orbit, fi is always negative; 
c total angular momentum; 

G z component of the angular momentum, G is positive or negative according as 
the motion is direct or retrograde; 

- t time of perigee passage; 

& argument of perigee; and 

fi right ascension of the ascending node. 

Exact expressions for three canonical constants a l9 a 2 , and a , denoted by Izsak as 
fi, c , and 6, respectively, are determined from initial values of the coordinates and their 
derivatives (Reference 4). Numerical values of these a’s are used to determine a certain 
set of orbital constants a Q , e 0 , and i Q (the initial values of the semimajor axis of the orbit, 
the eccentricity of the orbit, and the angle of inclination, respectively). These are used to 
find the p lf p 2 , A, and B (the perigee of the orbit, the apogee of the orbit ; and coefficients 
in the quartic polynomial F(p) — see Appendix B — respectively) necessary to factor 



F ifl) ~ ~ 2a t (p- p,) (p 2 “ p){p 2 + Ap + B) t ( 4 ) 

where p = a[ l - e) and p = a{ l + e) . Similarly, this same quartic designated by Izsak as 
P (/>), Equation 3, can be factored into a form which is equivalent to that of F(p), Equation 4. 
That is, 


P(p) “ “2h (P 2 -P)(p-Pi) [(/° ~ aK ) 2 + a2 ^ 2 ] » 


and we find that 


p 2 + Ap + B 


ip- a*) 2 + a 2 \ 2 . 
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The values for p, + p 2 , p l p i , A, and b are determined initially by solving the follow- 
ing system of nonlinear equations: 

p i + Pi ~ A = ~P a i 1 = 2a o • ] 


n + P, P 2 - (p, + P 2 ) A - c 2 - \ a * a~ 1 - c 2 + a 0 p Q 

(p, + p 3 ) B - p, p 2 A = = 2a 0 c 2 , 

PjP 2 b = ~T c2 ( a 2 2 " a 3 2 ) a r' = a oPo c2sin2i o • 


>■ 


J 


Vinti (Reference 4) has given a second-order solution of this system by a method of 
successive approximations. However, if higher order accuracy is desired, it is first 
necessary to obtain additional terms in the series solutions; this is a laborious task. A 
numerical method to obtain the solution is given in the next section. 


NEWTON-RAPHSON ITERATION SCHEME 

The solution of a set of nonlinear algebraic equations usually involves a great deal 
more work than that needed for linear systems. When n, the number of equations, is large, 
the solution of linear systems entails considerable computation time even on high-speed 
computers; the solution of nonlinear systems may often be almost prohibitive. 

The Newton- Raphson method (References 5, 6, and 7) can easily be applied when a 
solution is required for only a few equations. 

To solve a system of nonlinear equations such as 

(Pj + P 2 ) - A = 2a 0 . (5) 

B + P 2 P 2 ~ (P, + P 2 ) A = c 2 + a 0 p 0 • (6) 

(Pl + p 2 ) B - p, p 2 A = 2a 0 c 2 . ^ 

Pj p 2 B = a Q p Q c 2 sin 2 i Q . 


(8) 
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with unknowns (p t + p 2 ) , p x p 2 , A, and B by the Newton-Raphson iteration we begin with a 
trial vector 

X (k + 1) = X (k) _ j-1 F ( x oo) , 



Denoting Equations 5, 6, 7, and 8 by f,, f 2 , f 3 , and f 4 , respectively, that is, 


f l( x l- 

x 2 , 

X 3 1 

x a ) 

= (^1 + ^2)~ A_ 

2a 0 - 

f 2 ( X 1 > 

X 2 * 

X 3 ’ 

X 4 ) 

= B + p 1 p 2 - (P, 

+ p 2 ) A - c 2 - a 0 p 0 

f 3 ( X l' 

x 2> 

X 3 , 

%) 

= (p, +p 2 ) B - P, 

P 2 A “ 2a 0 c 2 , 

f 4 ( x l, 

x 2 . 

X 3 ’ 


= Pj P 2 B - a 0 p 0 

c 2 sin 2 i fl , 


we introduce the usual Jacobian matrix 



It will be noted that for the initial vector, the Jacobian determinant can be written 


10-1 0 
0 1 -2a 0 1 

0 0 -a 0 p 0 2a 0 

0 0 0 a 0 p 0 




The condition that S J i is not close to zero will be satisfied provided e 0 is not close to 
unity. 



8 


Next, we determine the exact inverse of Equation iO; only the final result is given 


here: 




11 Z 1 2 Z 13 Z 14) 


7 7 7 7 

21 22 23 ^24 


7 Z 7 Z 
31 ^32 33 ^34 


V 


7 7 7 7 

^4 1 ^4 2 ^4 3 ^4 4 / 


where 


Z n = 1 - 


(p, P 2 - B) (B - A 2 ) - BA (a +Pj +p 2 ) 

A 


Pj P 2 A + (Pi + P 2 )B 

g— 


Pi P 2 " B 


A + Pj + P 2 


'14 


Z 2. = A “ 


(pj P 2 ) (A+P, +p 2 ) (B-A 2 ) - BA(A+P t + P 2 ) 2 + Ba[(p! P 2 - B) + A(A + p t + P 2 


f (/°1 ^2 -b)[-a(a+p 1 + p 2 ) + b] - B(A+p, +p 2 ) : 

1 - 1 k A- 


Z 2 3 ~ (^1 ^2) Z 1 


(Pj P 2 )(A+Pj +p 2 ) 


■(a + p, + P 3 ) 2 + (p, P 2 -B) + A (A +P t + P 2 ) 


"2 4 


3 1 


^ii' 1 ) = 


'(Pj P 2 -B)(B-A S ) -BA(A+Pj+p 2 ) 

A 


7 =7 

32 1 2 


K P 2 )A + (Pj +p 2 )b 

2 


7 =7 

^33 ^13 


(pi p 2 " B ) 


Z 34 = Z 1 4 


(A+P! +P 2 ) 


b2 (Pi +P 2 ) + BA(P,P 2 ) 


Z 4 1 BZ 3 2 2T 
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Z 4 2 = BZ ,3 


B ('°l P 2 ~ B ) 

A 


Z 43 = BZ 34 


' B (A+Pj + P 2 ) 
A 


Z 44 = " Z 33 + ^34 


(p, Pj -B) + A (A + Pj +p 2 ) 

s 


where 


A ' " [(^r B ) 2 + A(p t p 2 -B)(A+p, +p 2 ) + b(A+Pj +p 2 )2] 


The Newton- Raphson iteration can now be written 




A solution will have been obtained when 


Max 


(x.< k + i> - XX k ) ) < 

i 

where e is any tolerance sufficiently small to obtain the 


e , 

degree of accuracy desired. 


The first-order solution through k Q is obtained in one iteration by beginning with the 
zero-order solution of Equation 9. It is expected that four iterations will be sufficient to 
obtain the solution through o(k 0 3 ). 

Table 1 gives the solution to the system as computed on the IBM 7090 using a single 
precision floating point Fortran program. The solution of Equations 5 through 8 was ob- 
tained accurately to seven significant digits in four iterations withe = kt 7 . 

The initial values of the unknowns, together with other necessary constants were com- 
puted from orbital data of Explorer XI (1961 v). 
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Table 1 

Numerical Results Obtained with the Initial Conditions 
p t + p 2 - 2a Q = 15.0588664, 

p x p 2 = a Q p Q = 56. 2660106 , 

A = 0.0, 

B = 0.0, 

c 2 = J 2 a E 2 = 0.044029034, 

sin i Q = 0.474484778, 

e = 1.0 x 10" 7 . 


Iteration 

Number 

Pi + Pi 

p \P i 

A 

B 

1 

15,0497347 

56.1626263 

-0.00913084351 

+0.00991251186 

2 

15.0497213 

56.1624885 

-0.00914439194 

+0.00993078329 

3 

15.0497213 

56.1624877 

-0.00914439194 

+0.00993078336 

4 

15.0497213 

56.1624877 

-0.00914439209 

+0.00993078336 


We immediately obtain the values of the elements a and e from 


a 


Pi + P 2 
2 


e 



DETERMINATION OF CANONICAL CONSTANTS FROM 
INITIAL CONDITIONS 

If the initial conditions (denoted by zero subscripts) t 0 , x Q , y 0 , z 0 , x 0 , y 0 , and z Q are 
known, we can determine a complete set of canonical constants a, e, S, _ t, ft, and w essen- 
tial to Izsak’s second- order solution. The canonical constants have the following meanings 

a semimajor axis of the orbit, 
e eccentricity of the orbit, 

S sine of the inclination of the orbit, 

- t in the limit as c 2 -0 of Keplerian motion - t is the time of perigee passage, 
ft in the limit as c 2 -0 of Keplerian motion ft is the right ascension of the ascending 
node, 

^ a constant of integration. 
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We first give the method of determining a, e, and s. The following expressions are 
computed: 


+ 



r o f o - x o *o + y 0 y 0 + z o *» 


Po r 


0 0 0 


r n + c 2 <x n i, 


pi + c2a 0 2 


7 _0 r o f o + P_o z o 

pi + c2a 0 2 


where cr g takes the sign of z 0 . Next, the following expressions are computed: 


x i " 2 v o 


PP Q 

Jfi + c 2 cr o 2 ) 


< 0 


a 3 _ x 0 yo _ y 0 x o ■ 


a 3 2 + (~ a 0 r 0 f 0 + p 0 *o) : 

1 - CT 0 2 


- 2a J c 2 a 0 2 , 


y 


2 



= cos 2 i 


o ’ 


o ’ 


f\ - y 2 


sin i 
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2 

P 


p o J 


P 1 P, 


P, + P 2 = 2P 0 x" 2 [l " K 0 x 2 y 2 - K q 2 x 2 y 2 (2x 2 - 3x 2 y 2 - 4 + 8y 2 ) + • ■ •] , 

= P 0 2 x' 2 [l + K 0 y 2 (x 2 - 4 ) - K 0 2 y 2 (l2x 2 - x 4 - 20x 2 y 2 - 16 + 32y 2 + x 4 y 2 ) + • • •] 


a 


T [Px + p 2 ) 


4 ^i Pji 

[p\ *p*Y ’ 


p! +p 2 


f x - ( 1_eJ ) 


Vo 


( sin »o) 


■TK„x 2 


V 


;4 y 2 


(7y 2 



( 12 ) 


where v 0 ~ sin I = S. Alternately tj 0 may be computed from 


a* - 2 c 2 
7,0-2 = 2(a 2 2 -a 3 2 ) 

We now have determined Izsak’s elements a, e, and S from initial conditions; these 
elements are accurate through o(k 0 2 ) and are used as input to the orbit computation 
procedure. 

We are now ready to determine the remaining canonical constants - t, fi, and We 
set 4> = E = 0 whenever e = 0. If e / o, we determine E from 


l +■ 


8a x c 2 


{ a 2 ~ a 3 2 ) 

r 


i* - 2a x c 


E 


tan* 


p 0 r o f o + cV o *0 


_ y-2fi | f[p - a*) 2 + a 2 A. 2 (a-p # )J 
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since 


cos E 


a 


-_P _ 0 

ae 


sin E 


P_0 r o *o + c2g 0 *o 
! V- 2fi j[p~ aw) 2 + a 2 Au 1 


where P 0 2 is given by Equation 11, c is given by Equation 2, and 


a 2 K 2 



(TT^yr (i-s 3 ) (i-ss^) 



(13) 




V 2 (l ~ S 2 ) 

1 - e 2 


4 (i-s 2 ) 

(l-e 2 )^ 


-4S 2 -e 2 ) 


+ • * * » 


(14) 


- 2h 


£L_ 

a{ !+/<)* 


(15) 


The angle 4> is completely determined by 

cos E — e t 

cos * = 1 - e, cos E ’ (16) 

j/ 1 - e , 2 sin E 

si r> $ = 1 - e, cos E ’ < 17 ) 

where e m is given by 

e - = e { 1 + 7~7( 1-2S2 ) + (TT^yi[(3-16S 2 + 14S-') - 2 (l ~ S 2 ) 2 e 2 ] + • • \ . (18) 

We next determine the angle '/'from 


0 , 

sin \p = “§“ ' S f 0 


(19) 


7 o r o f o + Po' z o 


i 2fi Vl - e 2 j /* 2 — - |/l - £ 2 sin 2 i/> 


cos \// 


( 20 ) 
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where (l - e 2 ) is given by Equation 12 and 

i 

(' - ' J ) 


«■ - 


1 - e 2 


k 2 + X 2 , 4^ 2 S 2 

: — = S 2 + 


(T^fpl 1 - 52 ) + -(T^V r(l " s2 ) [( 1_3S2 ) ■ ( 1 + sJ ) e2 ] + 

Next we compute a f, mean anomaly” M from 

M - E “ Kj e sin E - K 2 <£ - K 3 sin c£ + K 4 sin 2<£ + K 5 0 “ K 6 sin 20 + K 7 sin 40 


(21) 


(22) 


where 


= i + ^iL_s!) S 2 (3 + e2) 

. _ „2 (l -e 2 ) 3 ' ' 


1 - e 2 


K, = 


(23) 


’2U ~ V r S ’ - [(24 - 96S»<78S«) - (8-llS>) S>.>] . (24) 

'^iifr< 4 - ss ‘) s '’' ( 25 ) 


3v 


4 Vl - e 2 S 4 e 2 
32(l -e 2 ) 3 




5 [ 2 (l ~ e 2 ) 16 ( 


„2 - e 2 v 4 yi - 

4 (l - e 2 ) 8 ( 


^^[(24 - 27S 2 ) - (8-llS 2 )e 2 ]|s 2 
r^T[( 6 -7S 2 ) - (2-3S 2 )e 2 ]|s 2 , 


v 4 /l - e 2 

K? = 64 (l - e 2 ) 2 S • 


We can now compute - t as follows: 


(26) 


(27) 


(28) 


(29) 


- t = — - t . 
n 
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where 


n 



3W1-S 2 ) 

2 (l ~ e J ) 


■ 3 v 4 (l-S») 
8(l-e»)» 


[(l + 11S 2 ) 


(l ~5S 2 ) e 2 ] 


( 30 ) 


The right ascension df the satellite a is determined from 


V^o 2 +c2 i 1 “ CT 0 2 


y 0 



When the right ascension a is known, the right ascension of the ascending node fi is 
computed from 


fl 


a 


tan" 1 


(vi - s 2 



+ R, </. 


R 2 sin 2^ + R 3 4 > 


+ R 4 sin <£ + R s sin 2<£ " R fi sin 3^ “ R ? sin 4<£ , 


where 


_ Vi - s 2 _ v 4 Vi - s 2 

1 2 (l - e 2 ) 16 (l _ e 2 ) 3 


[(30 - 35S 2 ) + (2 + 3S 2 ) e 2 ] 


( 31 ) 


3^ 4 Vl - S 2 
* 2 ' 32(1 -e 2 ) 2 


S 2 , 


( 32 ) 


i]rr ^ (2 + e 2 ) 4 ^('-e 2 ^ [( 24 - 56SJ ) " ( 4+ MS 2 )e 2 - ( 2 + 3 S 2 ) e 4 ] ,(33) 


y* Vl - S 2 


2v 2 Vl - S 2 v 4 Vl - S 2 


(l - e 2 ) 2 4 ( 


rr^T [(4 ~ 28S 2 ) - ( 6 + 7S 2 )e 2 ]J e , 


( 34 ) 


V 2 Vl - S 2 i> 4 Vl - S 2 r , ,11 , 

4 (l - e 2 ) 2 " 8 (l “ e 2 ) 4 [“ + ( 1+S )* ]f 6 ’ 


( 35 ) 
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R 


6 


v 4 i i - s 2 

4(l -e 2 ) 4 


(2 -s’) 


( 36 ) 


R 


7 


v 4 i 1 - s 2 

64 (l ~ e 2 ) 4 


(2 + S 2 ) 


( 37 ) 


We next compute w and v, which are analogous to the argument of latitude in Keplerian 
motion and the true anomaly in Keplerian motion, respectively, from 

W = 4* - Mj sin 2^ + 3M 2 sin 4</^ 


and 


V - <p + L t sin 2<£ + L 2 sin 4 , 


where 


( 38 ) 


M 


1 


8 



e 4 

M 2 “ 256 ’ 



3k 4 

2 “ 256 ’ 


The mean argument of perigee co is given by 

co - w - V . 


The constant of integration co is given by 

co - W - (1 + e)v 


where 


€ 


2 4 

4( 1 ! e ^)^ 12 " 15S2) + 64 {y~ e 2 ) * S 288 " 1296SI + 1035S4 ) 
- (l44 + 288S 2 - 510S*) e 2 ] + ••• . 


( 39 ) 
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ORBIT COMPUTATION PROCEDURE 

With the exception of the velocity formulation, the computational procedure developed 
here makes use of the unmodified expressions of Izsak (Reference 1). 


Input 

The 11 inputs are: a, e, s, - t, A, a;, j 2 , a E , At, and T f . The six inputs a, e, s, - t, 
fi, and ^ are constants of integration (see definitions on page 10). The other Inputs have 
the following meanings: 

J 2 the coefficient of the second-degree Legendre polynomial in the earth’s gravita- 
tional potential, 

a E the earth’s equatorial radius, 

ji the product CM where G is the gravitational constant and m is the earth’s mass, 

At time interval of integration, 

T f final time. 

The following values of j 2 , and a £ determined by W. M. Kaula (Reference 8) were 
used in the computations: 

h- = 3.986032 x 10 2 megameters 3 ksec" 2 , 

J 2 = 1.0823 x lO' 3 

a E - 6.378165 megameters. 


Equations and Fundamental Constants 

From Vinti's expression (Equation 2) and the input constants determined by Kaula, we 
have c = 0.20983097 megameters. In addition to a 2 \ 2 (Equation 13), * (Equation 14), - 2h 
(Equation 15), costf (Equation 16), sin<£ (Equation 17), e. (Equation 18), t 2 (Equation 21), 
(k 2 +\ 2 )/v 2 (Equation 22), n (Equation 30), V (Equation 38), and e (Equation 39), the follow- 
ing equations are used in the computation: 


k 2 = 


1/2 g 2 j;4 g 2 p -» 

52 - - 1052 + 11S< ) + s< e2 J + 
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and the generalized Kepler equation 

E - Kj e sin E = n (t - f) + K 3 <f> + K 3 sin<£ “ K 4 sin 2<f> ~ K 5 0 + sin 2^ - K ? sin # * (40) 

where the K. are given by Equations 23 through 29. 

The right ascension a is computed from 

a = Cl + tan 1 (Vl - S 2 tan - F + F 2 sin “ R 3 <£ “ R 4 sin 

“ R s sin 2<j> + R 6 sin 3$ + R ? sin 4<£ t 

where the R i are given by Equations 31 through 37. 

The argument of latitude \p is computed from the following equations: 

W = (l + e)v + oj , 

4> = W + Mj sin2W + M 2 sin4W + ••• . 

The mean argument of perigee co is computed from 

co = eV + co . 

The anomalistic mean motion is computed from 


* * 


The motion of the node 77 is computed from 


V 


3v 2 Vl - S 2 

2 (l - e^) 2 


3v« Yl - s 2 

16 (l ~e J )« 


[(l8~ 13S 2 ) 


+ 24S 2 e 2 ] - • • ■ . 


The oblate spheroidal coordinates is computed from 

a - s sin \p . 
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The z component of the angular momentum G is computed from 


G 2 


-2fla 2 (l -S 2 ) 




The oblate-spheroidal coordinate p is computed from 

p = a (l - e cos E) . 

Initially for t = t 0 we start with values <£ and determined from Equations 16, 17, 19, 
and 20 to solve the generalized Kepler equation given by Equation 40 using a Newton- 
Raphson iteration scheme. We test |e(<£. + 1 , ^ i + 1 ) -e(<£., ^.) | <€, where e > o was chosen 
to be 10' 7 . In general, only two or three iterations are required before sufficiently ac- 
curate values of e, <£, and are obtained. The oblate -spheroidal coordinates p, a, and a 
are then computed; p, <?, and a are then used to calculate x, y, and z. 


OUTPUT 

This program generates position and velocity for equally spaced intervals of time. 
Oblate- spheroidal coordinates are defined by the equations 

x “ ip 2 + c 2 ii ~ a 2 cos a « 
y = ip 2 + c 2 / 1 “ a 2 sin a * 


z - pa 


t ~ y p 2 + c 2 [l -cr 2 } . 


The formulas for velocity are given in Appendix B 5 they are 


x - - ay + x 


pp aa 


P 2 + c 2 1 - a 2 


y = + ax + y 


/ PP Qg \ , 

+ C 2 1 “ cr 2 y 


z - pa + ap 



20 


where 


P - 


— ~ ^ - ae + a 2 \ 2 )sinE * 


p 2 + C z <7 




p 2 + c 2 cr 2 


a Vl - e 2 \ Vl - t 2 sin 2 ^ cos ^ 

r 


(p 2 + c 2 ) (l - a 2 ) 


REMARKS 

The computational procedure as it exists in this report was programmed by the author 
in single-precision floating-point Fortran for an IBM 7090 computer at the Goddard Space 
Flight Center, All machine results were compared with hand calculations and the practi- 
cality of the method was confirmed. The procedure is presently being compared with both 
single and double precision numerical integration. 
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Appendix A 

List of Symbols 


A 

a 



B 

c 


c 


E 

*(*,• *.) 


E (^i+i’ ^i+0 


e* 

F 


Fip) 

f 2» f 4 


coefficient in the quartic polynomial Tip). See Appendix B. 

canonical constant, one of Izsak's elements, semimajor axis of the orbit. 

the earth's equatorial radius. 

initial value of the canonical constant a. 

coefficient in the quartic polynomial F(/o). See Appendix B. 

a constant defined by Vinti's expression c 2 = j 2 a E 2 . 

a canonical constant; in the limit as c 2 -0 of Keplerian motion c is the 
total angular momentum. 

angle corresponding to the eccentric anomaly. 

the ith value of the eccentric anomaly. 

the ( i + 1) th value of the eccentric anomaly. 

canonical constant, one of Izsak' s elements, eccentricity of the orbit. 

initial value of the canonical constant e. 

second eccentricity. 

the earth’s force function. 

quartic polynomial fundamental to Vinti's theory. 

representation of a set of four equations to be solved by the Newton- 
Raphson method. 


G the gravitational constant. 

G a canonical constant; in the limit as c 2 -0 of Keplerian motion G is the 

z component of the angular momentum. G is positive or negative ac- 
cordingly as the motion is direct or retrograde. 



the Hamiltonian. 


q n ; Pj , p 2 , ■ ■ • , p n ; t) the Hamiltonian function (in which time appears 

explicitly) of a dynamical system of n degrees of freedom with n gen- 
eralized coordinates , q 2 f • ■ • , q n and the conjugate momenta 

Pi* P 2’ P n * 

a canonical constant; in the limit as c 2 - 0 of Keplerian motion h is the 
total energy in the orbit and always negative. 

one of Izsak T s elements, inclination of the orbit. 

angle of inclination. 

initial angle of inclination. 

the Jacobian matrix of the Newton-Raphson method, 
the Jacobian determinant. 

the coefficient of the second-degree Legendre polynomial in the earth’s 
gravitational potential. 

notation used for the coefficients in Kepler equation. 

the value . 

p o 2 

modulus appearing in elliptic integral of the first kind, 
modulus appearing in elliptic integral of the first kind, 
the Lagrangian. 
the earth’s mass. 

"mean anomaly". 

the anomalistic mean motion. 

a constant used in the generalized Kepler equation, the auxiliary mean 
motion. 
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Po 

P(/0) 

P/3 1 Per ' Pa 

Q(cr) 

R i 


r 


S 


T f 

t 

*0 

At 

- t 


U 

V 

V 
v 


w 

w 

X 


x, y, z 
x 0’ yo’ z o 


the value . 

quartie polynomial fundamental to Vinti’s theory, 
the generalized momenta. 

quartie polynomial fundamental to Vinti’s theory. 

notation used for the coefficients in the equation for right ascension of the 
ascending node. 

the geocentric distance of the satellite, 
the initial geocentric distance of the satellite. 

canonical constant, one of Izsak’s elements, sine of the inclination of the 
orbit. 

final time. 

time. 

initial time. 

time interval of integration. 

a canonical constant; in the limit as c 2 - 0 of Kepler ian motion - t is the 
time of perigee passage. 

the speed of the satellite. 

a ’’true anomaly” analogous to that in Keplerian motion, 
the potential which Vinti obtained in oblate spheroidal coordinates, 
velocity of the satellite, 
initial velocity of the satellite. 

’’argument of latitude” analogous to that in Keplerian motion. 

Hamilton’s characteristic function. 

a trial vector for the solution of a set of nonlinear equations by the 
Newton- Raphs on method. 

coordinates in the rectangular system. 

the initial values of the coordinates in the rectangular system. 
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X, y, z 

*o- *0- i o 

z ii 

a 

V a 2 , and * 3 
§ 


the velocity coordinates in the rectangular system, 

the initial value of the velocity coordinates in the rectangular system, 
element of the inverse Jacobian matrix, 
the right ascension of the satellite. 

Vinti's canonical constants denoted by Izsak as h, c, and G 
respectively. 

the declination of the satellite, 
the motion of perigee. 


€ 

v 

7} 0 % sin I = s 


K 






P, v, a 

°V a o 


p , CT, a 


4 > 

* 

fl 


an arbitrarily chosen small positive real number (used as a tolerance 
in the Newton- Raphson method). 

the motion of the node. 

a series used in the computation: defined by Equation 14, 

the product GM where G is the gravitational constant and M is the earth’s 
mass. 

a dimensionless parameter of the order 10" 3 in the case of the earth. 

coordinates in the oblate spheroidal system. 

the initial condition of the coordinates in the oblate spheroidal 
system. 

the velocity coordinates in the oblate spheroidal system. 

the initial conditions of the velocity coordinates in the oblate 
spheroidal system. 

perigee of the orbit. 

apogee of the orbit. 

"true anomaly". 

"argument of latitude". 

a canonical constant; in the limit as c 2 - o of Kepler ian motion fi is the 
right ascension of the ascending node. 


a canonical constant, one of Izsak’s constants, a constant of integration. 

a canonical constant; in the limit as c 2 - Oof Keplerian motion £ is the 
argument of perigee, 

the mean argument of perigee. 
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Appendix B 

Derivation of the Velocities 
in Vinti’s Accurate Intermediary Orbit 
of an Artificial Satellite 


Introduction 

Izsak (Reference Bl) has given an analytic solution for Vinti’s intermediary orbit, with 
both periodic and secular terms correct through the second order in a certain oblateness 
parameter y = c/a (to be defined later). His solution giving the position vector of the satel- 
lite makes extensive use of Jacobian elliptic functions, linear transformations, and map- 
pings in the complex plane. Vinti (Reference B2) also has given an analytic solution to this 
problem of satellite motion using rapidly converging infinite series instead of Jacobian 
elliptic functions. His solution not only gives the periodic terms correct to the second 
order, but also the secular terms to an arbitrarily high order. 

This appendix presents the derivation of the velocity vector through the use of equa- 
tions from both Vinti and Izsak. However, the orbital elements used in this derivation were 
introduced by Izsak. 


Determination of Velocity 

The oblate spheroidal coordinates p , <j, and a are defined by 


ip 2 + c 2 

H - or^ cos a , 

(Bl) 

ifp 2 + c 2 

i 1 - cr 2 sin a , 

(B2) 

pa , 


(B3) 


r 


i + c 2 (1 - <t 2 ) , 


(B4) 
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where a is the right ascension of a satellite; r is the geocentric distance; and c is a con- 
stant defined by c 2 = J 2 a E 2 . The quantity J 2 is the coefficient of the second-degree 
Legendre polynomial in the earth's gravitational potential 

00 

V = T 1 P n< sin8) 

*- n = 1 

(B5) 

where 8 is the declination of the satellite, a E is the equatorial radius of 
fi = GM, where G is the gravitational constant and M the mass of the earth. 

the earth, and 

Differentiating Equations B1-B3 with respect to time we find 


• _ . . ( PP CTCT \ 

x - -«y + *^ 2 — +c2 • 

(B6) 

• _ . . , ( PP <T& \ 

y ' +ax + y [p> + c 2 ■ 1 - W ’ 

(B7) 

Z = p <J + <jp . 

(B8) 


Squaring and adding Equations B6-B8 we obtain 


u 2 


y 2 


+ z * 


/ p 2 + c 2 a 2 \ 
\ p 2 + c 2 / 


P 2 


+ 



+ Q l (J 




2 + 


{P 2 


+ C 2 )(l 




(B9) 


The expressions for p, a, and a can be obtained from the following equations, which define 
the generalized momenta: 


= 


dh 

dp 

. , . _ dS _ . VT(p7 

h J P dp p2 q2 

(BIO) 

dL 

do- 

t 3 . dS Vo(o-) 

h 2 a " da - ± J _ a 2 ' 

(Bll) 

dL 

da 

, , . as 

h 3 a " da - G ’ 

(B12) 


*The caret above h, c, and G (that is, fi, c, and G) indicates canonical constants, referred to by Izsak (Reference Bl), where 
£ is the total energy in the orbit and always negative, c is the total angular momentum, and 6 is the z component of the 
angular momentum, posicive for direct motion. 
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where 




p 2 + C 2 Cr 2 

p 2 + C 2 


h / 


p 2 + C 2 a 2 ^ 2 

1 - a 2 ’ 3 


(p 2 + c 2 ) (l - a 2 ) * 


P(p) = 2hp 4 + 2pp 3 - (c 2 - 2c 2 h)p 2 + 2 c 2 pp - c 2 (c 2 ~ G 2 ) , 
Q(cr) = - 2c 2 h a 4 - (c 2 - 2c 2 h)a 2 + ( c 2 - G 2 ) - 


(B13) 

(B14) 

(B15) 


s = Sip, er, a) is the action function, and l is the Lagrangian given by L - T - V, where 


T = 


1 /ds\ 2 

2 \dt/ 


■^^hj 2 p 2 + h 2 2 a 2 + h 3 2 a*j 


V = V(p, cr) = 


- PP 

0 2 + C 2 cr 2 


Here ds/dt is the speed along the path and V(p, <?) is the potential function introduced 
by Vinti (Reference B3), 


The radicand in Equation BIO can be written in the form 


Pip) = - 2h (p - p x ) (p 2 - p)(p - p 3 ) [p - p 4 ) , (B16) 

where p 1? p 2 , p 3 , and p 4 are the zeros of P(p). Izsak (Reference Bl) has given the zeros 
in the form 


P x - a(l - e), p 2 - a(l + e], p 3 = a {k - iX), p 4 - a [k + iX). (B17) 

The orbital elements a and e are the semimajor axis and the eccentricity of the orbit, re- 
spectively, even though it is not an exact ellipse. They are defined by the first two of Equa- 
tions B17. The quantity i is the imaginary unit (/-*!). 

If we substitute p 3 and p 4 from Equations B17 into Equation B16: 

Pip) = -2h(p - Pj) (p 2 - p) [p 2 - 2a*p + a 2 (* 2 + X 2 )] . 

The quantities * and k 2 +X 2 are given in Reference Bl in terms of a, e, and s = sin i, 
where I is the inclination of the orbit: 

y 2 (l ~ s 2 ) (l - e 2 ~ 7 2 s 2 ) 

(l - e 2 - y 2 ) (l - e 2 - y 2 s 2 ) + 4y 2 s 2 


K 


(B18) 
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s . x2 y 2 s 2 [(l -’e 2 - y 2 ) (l - e» - y 2 s 2 ) + 4y 2 
(l _ e 2 - y 2 ) (l _ e 2_ y 2 s 2 ) + 4y 2 s 2 

where y = c/a, a small dimensionless parameter. 

The quartic Q[a) contains only even powers of o- and can be written 


(B19) 


Q(cr) = _ 2c 2 h i ~ °" 2 ] (&2 “ a2 ) * 

where the four real zeros of Q(o-) = 0 are to^and ±a 2 : 0 < & l < 1, <x 2 » l. As pointed 

out by Izsak (Reference Bl), o- oscillates between the values -o*, and + o- x . Therefore, ct-j 
is a convenient parameter to use as the sine of the inclination I of the orbit. 

When we introduce Izsak f s formulas, 


p x - a ( 1 - e), p 2 - all + el , 
p - a ( 1 - e cos E) , 


<j 1 - s sin I 


a - s sin \p , 



c 4 s 2 — - a 4 (l - e 2 )(* 2 + \ 2 ) , 


where E is the eccentric anomaly and \p is the argument of latitude, and several of the afore- 
mentioned relations into Equations BIO and Bll we obtain 


P 


p 


V- 2h 

p 2 + c 2 


ae 


/ p 7 


2a Kp + a 2 Jp + X 2 ] sin E 


(B20) 


i - 2h 


/l - e 2 }^k 2 + X. 2 ]/l - Z 2 sin 2 i/i cos i// . 


1 - CT : 


(B21) 
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The coefficients and the roots of the quartic equation p(p) =0 can be related to those 
of Q(cr) = o in the following manner: 


2a ( 1 + k) 



(B22) 


a 2 1 - e 2 ) + 4* + (k 2 + X 2 )J 


- 2c 2 h 
2h 



2a 3 - e 2 ) k + (k 2 + \ 2 )J “ ” c 2 ~^r f 

a«(l - e 2 ) (k 2 + \ 2 ) = -c 2 ~ ~ ° 

2h 


s' 

u 


(B23) 


(B24) 


(B25) 


Consider the following expression for l 2 given by Izsak (Reference Bl) 

(l - e 2 - y 2 ) (l - e 2 - y 2 s 2 ) + 4y 2 s 2 
(l - e 2 “ y 2 ) (l - e 2 - y 2 s 2 ) + 4y 2 


y 2 s 2 


(B26) 


If we substitute c/a for y and solve for s 2 /l 2 we obtain 

s 2 _ a 2 (l - e 2 ) f [a 2 (l - e 2 ) - c 2 ] [a 2 (l ~ e 2 ) - c 2 s 2 ] + 4a 2 c 2 \ 
l 2 " c 2 [[a 2 (i - e 2 ) ~“c*] [a 2 (l - e 2 ) - c 2 s 2 ] + 4a 2 c 2 s 2 j 


(B27) 


Next we introduce a parameter p = a(l -e 2 ), the semilatus rectum, which Vinti (Refer- 
ence B2) uses in his oblateness parameter k = c 2 /p 2 . It is clear that Equation B27 can be 
written 



ap (ap ~ c 2 ) (ap * c 2 s 2 ) + 4a 2 c 2 
c 2 (ap - c 2 ) (ap - c 2 s 2 ) + 4a 2 c 2 s 2 


(B28) 


From Equation 4.13 of Reference B2: 

_ 2 2 (ap - c 2 ) (ap - c 2 V) + 4a 2 c 2 

0 (ap - c 2 ) (ap - c 2 7j 0 2 ) + 4a 2 c 2 Tj 0 2 


(B29) 



34 

where v 0 = s = sin i, we see that c A s 4 // 2 = Bap- Solving for c 2 s 2 /l 2 from Equation B23 and 
inserting it into Equation B25, we obtain 

G 2 - (l - s 2 )(c 2 + 2c 2 hs 2 ) (B30) 

From Equation B23, 

— = -c 2 (l - S 2 ) + . (B31) 

2h 1 

Multiplying Equation B28 by c 2 to obtain c 2 s 2 /l 2 and inserting it into Equation B31 we find 

(B32) 

Vinti has given an identical result in Equation 4,15 of Reference B2: 


2h 


- -c 2 (l~s 2 )+ap 


(ap - c 2 ) (ap - c 2 s 2 ) + 4a 2 c 


(ap - c 2 ) (ap - c 2 s 2 } + 4a 2 c 2 s 2 J 


Tsf - s a «Po = _c2 i 1 - V) + a P 

where a 0 is a semimajor axis, p 0 the semilatus rectum. Since h = a lt 6 = a 2 , and G = a 3 , 
we can easily rewrite Equation B30, using Equation B33, to obtain the final result for 6, 


(ap 


2 ) (a p - c 2 V) 


+ 4a 2 c 


2^2 


(ap “ c 2 ) (ap - c 2 Vq) + 4a 2 c 2 ?7 0 2 J 


, (B33) 


a nPr 


(1 - .») . 


(B34) 


Equation B34 is equivalent to Vinti’s Equation 4.15a of Reference B2, 


a 3 = a 2 [1 - 


c 2 t? 0 2 \ 2 


a oPo 


cos I. 


Using Equation B34 we obtain, as in Equation B12, 


G. 


(B35) 


It should be noted that the following formulas relate k and * 2 + k 2 to Vinti T s A andB: 

a 2 (k 2 + k 2 ) = B , (B36) 


- 2a k 


A 


(B37) 
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where A and B are given by Equations 4.12 and 4.13 of Reference B2: 


~2ac 2 (l - rj 2 ) (ap - c 2 r ?Q 2 ) 
(ap - c 2 ) (ap - c 2 77 q 2 ) + 4a 2 c 2 77 0 2 


B ~ C 2 T) 2 


(ap - c 2 ) (ap - c2 Vq) + 4a 2 c 2 
0 (ap - c 2 ) (ap - c 2 ?j 0 2 j + 4a 2 c 2 77 0 : 


From Vinti’s Equation 4.16 (Reference B2) we have 


c 2 (ap - c 2 ) (ap - c 2 t? 0 2 ) + 4a 2 c 2 7? Q 2 

a P (ap “ c 2 ) (ap - c 2 r] 2 ') + 4a 2 c 2 


c 4 77 2 (Bap) 


where t ? 2 = ct 2 . Using Equations B36-B38, together with the values for P (p) and Q(o-), we 
find that in Vinti’s notation the quartics P(p) and Q(a) can be factored in the form 

Ftp) “ - 2a x (p 2 - p) (p - Pj) (p 2 + A p + B) , (B40) 

G(t?) = ~ 2a x c 2 (77 0 2 - 7) 2 ) (?) 2 2 - 7? 2 ), (B41) 


The following equations for p, 6 ■, and a are easily obtained from Equations B10-B13, 
B20, B21, and B35: 


p = — “ r ae ip 2 ~ 2a kd + a 2 (K 2 + X 2 ) sinE 

p 1 + c ' ' 


— ~ ~r / 1 - e 2 V k 2 + k 2 V 1 - l 2 sin 2 \p cos 1 p ^ 

p l + c i cr 1 ^ 


(p 2 + c 2 ) (l - a 2 ) 


If we write the equations for p, <7, and a in Vinti’s notation, we obtain 


P - — ; ae Vp 2 + A p + B sinE, 
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V 


V - 2ttj 

p 2 + c 2 tj 2 


CV 0 V 2 


VT 


q 2 sin 2 i p cos \p. 


(B46) 


4> 


(P 2 + C 2 ) (l - 7? 2 ) 


(B47) 


where q = i 7 0 /tj 2 . 

Now, substituting Equations B42-B44 into Equations B6-B8 we obtain the following: 
. _ - Gy xa V' 1 - 2h 

( p 2 + C 2 ) (l - c r 2 ) p 2 + C 2 Cr 2 


pe jp 2 - 2axp + a 2 (* 2 + k 2 ~) 


+ c J 


sin E 


a 


/T - e 2 
1 - a 2 


a 

c 


V^TT 2 


iT 


l 2 sin cos s/j 


(B48) 


y - 


+ Gx 


ya 


VTaT 


[p 2 + C 2 ) (l - <T 2 ) p 2 + C 2 cr 2 


pe 


ip 2 - 2a*p + a 2 (x 2 + \ 2 ) 


d 2 + c 2 


sin E 


a V 1 ~ e 2 

1 - cr 2 


a 

~c 


i k 2 + k 2 y 1 - l 2 sin 2 i// cos 0 


, (B49) 


a V- 2h~ 

P 2 + C 2 £T 2 


~ p y 1 - e 2 y x 2 + \ 2 V 1 “ Z 2 sin 2 \p cos i/? 


+ ere V~p 2 - 2a*p + a 2 (k 2 + \ 2 ) sinE 


(B50) 


The velocity components given above are now being used in an orbit determination 
program formulated by the author* 
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